arXiv:1508.05036v2 [cond-mat.supr-con] 25 May 2016 


Population imbalanced lattice fermions near the BCS-BEC crossover: 
II. The FFLO regime and its thermal signatures 


Madhupama Karmakar and Pinaki Majumdar 
Harish-Chandra Research Institute, Chhatnag Road, Jhunsi, Allahabad 211 019, India 

(Dated: May 26, 2016) 

We study the effect of high population imbalance in the two dimensional attractive Hubbard model, in the cou¬ 
pling regime corresponding to BCS-BEC crossover, and focus on thermal effects on the Fulde-Ferrell-Larkin- 
Ovchinnikov (FFLO) state. Using a method that retains all the classical thermal fluctuations on the FFLO state 
we estimate a very low Tc and infer a strongly first order normal state to FFLO transition. The Tc is an order of 
magnitude below the mean field estimate. We track the fermionic momentum distribution, the density of states, 
and the pairing structure factor deep into the normal state. The pairing structure factor retains weak signature of 
finite momentum pairing to a high temperature despite the low Tc itself, while the spin resolved density of states 
changes from the ‘pseudogapped’ FFLO character to gapless and pseudogapped again with increasing temper¬ 
ature. These results map out the rather narrow temperature window, and diverse physical indicators, relevant to 
FFLO states on a cold atom optical lattice, and complements our work on the lower field ‘breached pair’ state. 


I. INTRODUCTION 

The last decade has seen an intense search for spatially 
modulated superconducting (or superfluid) states which were 
originally predicted by Fulde and Ferrell (FFj^ and Larkin 
and Ovchinnikov (LOj^. These modulated paired states are 
expected to arise in the presence of a strong magnetic fleld 
that Zeeman couples to the fermions, but in most traditional 
superconductors the ‘orbital’ suppression, via Abrikosov lat¬ 
tice formation, takes place before any spin related effects can 
show up. This led to a decades long wait before a heavy 
fermion, 13221^ CeCoIns, som e Qua si-2D organic superconduc¬ 
tors of the /^-BEDT famil}!^^®, and the ferropnictidJ^®^ 
KFe 2 As 2 were identifled as promising ‘solid state’ candidates. 
Developments in cold atom physics, on the other hand, al¬ 
lowed the exploration of pairing among population imbal¬ 
anced neutral fermions, without the complication of orbital 
(Lorentz force) effectJ^^®^. 

The modulated, or flnite momentum, character of the paired 
state is difficult to probe since external flelds do not cou¬ 
ple directly to the pairing order parameter. Indirect probes 
in the solid state, for example, have probed (i) t herniody- 
namic features like speciflc heat and magnetizatioiP^I^M^ 
(ii) spectroscopic aspects like NMR shift and relaxatioiP®^, 
or (iii) spatial modulation of the induced magnetization, and 
have yielded suggestive results for the high fleld state in 
CeCoInJ^5C2l^ the organicJ^^®^, and the pnictideJ^^BUl. in cold 
atomic gases the direct spatial signature of a modulated super¬ 
fluid state continues to be elusive but the density profiles of the 
up and down spin condensates give clear indication of popula¬ 
tion imbalance^SED finite spin polarization indicateP^®!! 
the coexistence of unpaired fermions with the superfluid con¬ 
densate down to zero temperature. 

There is a large body of theoretical work exploring modu¬ 
lated mean fleld states in both the continuum and lattice cases. 
The mean fleld theory (MET) is non trivial since the nature 
of modulation, i.e, the wave vectors inv olved, is sensitive to 
the applied fleld or population imbalanc j^^ ^i Nevertheless, 
the density and fleld (or magnetization) window over which a 
FFLO ground state is expected is reasonably mapped ouPSEl 


in lattice models. In fact in two dimensions there now exists 
a detailed characterization of the FFLO ground state based on 
extensive variational calculation^^. 

What most of these calculations do not address is the ther¬ 
mal window over which the FFLO state is expected to sur vive. 
Mean fleld theory has been extended to flnite temperaturJ^^ESI 
but in the BCS to BEC crossover regime, where one expects 
the Tc to be highest, the results are very unreliable. In three 
dimensions, for couplings chosen to obtain the peak zero fleld 
Tc in the BCS-BEC crossover window, the EELO Tc is es¬ 
timated to be about half the zero fleld Tc, which mean fleld 
theory sets as ^ l.lt (where t is the hopping amplitude). This 
is a severe overestimatJSli Dynamical mean fleld theor}P2HlIl 
(DMET) has been employed, both in the ‘real space’ as well 
as the cluster form, but typically addressing only strongly 
anisotropic lattices. Einally, there is quantum Monte Carlo 
(QMC) data in two dimensionwhich sets the zero fleld Tc 
as ^ O.lt (and the EELO Tc obviously less than that) but sug¬ 
gests that EELO correlations survive to a scale T ^ Tc. 

As the numbers above indicate, estimates of the thermal 
stability window vary over a wide range. The existing results 
also do not address the spectral features of the thermal EELO 
state, or the normal state with EELO correlations. We feel 
there is room for a less computation intensive, possibly more 
intuitive, approximation that retains the rich detail of the mean 
fleld ground state but also accesses the key thermal fluctua¬ 
tions. 

We address the thermal physics in the EELO window, with 
the coupling set to the BCS to BEC crossover regime, in the 
two dimensional attractive Hubbard model. Using an auxil¬ 
iary fleld based real space Monte Carlo we observe the fol¬ 
lowing: (i) The maximum Tc for EELO phases is about 1/6 
of the zero fleld Tc, and about 20% of what mean fleld theory 
predicts. In absolute terms it is just about two percent of the 
hopping scale, (ii) The EELO to normal transition is strongly 
first order and the pairing structure factor retains a weak signa¬ 
ture of flnite momentum pairing to several times the Tc scale, 
(iii) The spin resolved density of states evolve from a ‘pseudo¬ 
gapped’ form, characteristic of LO order, at low temperature 
to a gapless form above Tc, but then loses weight at the Eermi 
level as the temperature is further increased. 











The paper is organized as follows. In Section II we briefly 
touch upon our model and method since most of this is dis¬ 
cussed in detail in Paper Section III presents our results on 
the variational ground state, and Monte Carlo inferred thermo¬ 
dynamic, spatial, and spectral features. Section IV discusses 
some issues and benchmarks related to our numerical method. 
We conclude in Section V. Much of the background material 
and the low held results are discussed in our Paper I. 

II. MODEL AND METHOD 

A. Model 

We study the attractive two dimensional Hubbard model 
(A2DHM) on a square lattice in the presence of a Zeeman 
held: 

H = Ho - - \U\'^ni-fnii (1) 

i i 

with. Ho = where Uj = -t only 

for nearest neighbor hopping and is zero otherwise. cFiz = 
(1/2) — n^). We will set t = 1 as the reference energy 

scale, /i is the chemical potential and h is the applied mag¬ 
netic held in the z direction, f/ > 0 is the strength of on site 
attraction. We will use Ujt = A. 

B. Methods 

Methodological issues have been discussed in detail in ear¬ 
lier papers^^EIl so we put in only a brief description for com¬ 
pleteness. The interaction is decoupled using a space-time 
varying complex auxiliary held A^(r) in the pairing chan¬ 
nel and we retain only the zero frequency mode of the aux¬ 
iliary held, i.e make a static auxiliary held (SAP) approxima¬ 
tion. However, the spatial fluctuations of are completely 
retained. The SAP scheme leads to the effective Hamiltonian: 

Heff = Ho-h'^ai^ + + h.c) + (2) 

i i 

where Hd = the stiffness cost associated with 

the auxiliary flelcP^. We use two strategies to study the model 
above: 

I. Monte Carlo: We generate the equilibrium {A^} conflg- 
urations by iteratively diagonalising the electron Hamiltonian 
Heff for every attempted update of the auxiliary flelds. The 
Monte Carlo is implemented using a cluster approximation, 
in which instead of diagonalising the entire L x L lattice for 
each local update of the A^ a smaller clust er, of size Lc x Lc, 
surrounding the update site is diagonalisecP^^. 

II. Variational calculation: The zero temperature limit 
within the SAP scheme is equivalent to unrestricted minimiza¬ 
tion of the ground state energy over conflgurations of the held 
A^. We have carried out minimization of the energy at several 
values of /i and h, exploring the following kind of periodic 


conflgurations: (i) ‘axial stripes’: A^ ^ Aq cos(gx^), diag¬ 
onal stripes A^ ^ Aocos(g(xi + ^i)), (ii) two dimensional 
modulations, A^ ^ Ao(cos(gx^) + cos{qyi)), and of course 
(hi) the unpolarised superfluid (USP) state A^ = Aq. We 
minimize the energy with respect to the g, and Aq. 

III. Green's functions for T = 0.* It is useful to set up a 
low order approximation for the Green’s function of the elec¬ 
tron, applicable in the T = 0 variational state, for Aq << zt, 
where the coordination number 2 : = 4 in 2D. In PPLO the Aq 
is strongly suppressed due to the magnetic held and the elec¬ 
tron pairing takes place between the k and — k±Q states. The 
resulting spin resolved Green’s functions take the approximate 
form: 

^-——- \ ^ . —r 

lUJn - (e(k) - jJiO - Stt(k, lUJn) 

/\ ^ 'Y 

I]n(k,*a;„) = x[(zu;„ + e(-k-Q)-Ai^) 

_ 1 _ 

(iLOn + e(-k + Q) - jJiO 

where, e(k) = —2t{cos{kx) -Vcos{ky)). A similar expression 
can be obtained for G^|(k, iujn) as well. 

One can extract the spectral function A||(k, cc) = 
— (l/7r)/mG^l(k, cc -f ip) \y^o from the above equation, 
whose k-sum would give the expression for the DOS. 


C. Parameter regime and indicators 

The results in this paper are at G = 4t, both within 
Monte Carlo and the variation scheme. We have also explored 
U = 2t variationally and observed that the PPLO regime is 
signiflcantly shrunk for a lower value of U/t. In fact, in order 
to access the PPLO regime at G = 2t we had to use a larger 
system size of L = 36 since at L = 24 (at which U = At cal¬ 
culations are being carried out) no PPLO signatures could be 
observed. Por the ground state, determined variationally, we 
have checked that the character of the /i — phase diagram 
does not vary much between L = 24 and L = 48. At 1/ = At 
we have explored the h — T dependence at multiple values 
of /i below half-fllling (the physics above half-fllling can be 
inferred from this) but the qualitative physics seems similar, 
so this paper focuses on detailed indicators at a single p. The 
density at this point is n ^ 0.94, and does not vary signifl¬ 
cantly for the horT that we have chosen. We discuss temper¬ 
ature dependence in the PPLO regime: h/t ^ [0.85 : 1.25]. 

As discussed in Paper I, the phases can be classifled into 
(a) unpolarized superfluid (USP), (b) breached pair (BP), 
(c) modulated superfluid (PPLO), and (d) a partially polarized 
Permi liquid (PPPL). 

We use the following indicators to characterize the physics 
on the L X L system: (i) The structure factors, Sa{^) and 
5'm(q) deflned earliei^S, and r(q) = 

where P^j = (ii) The bulk magnetization, 

and the pairing order parameter 5 'a (Q), where Q is the order¬ 
ing wave vector, (iii) Monte Carlo snapshots of (a) the magni¬ 
tude \Ai \ of the pairing held, (b) the correlation cos{0o — Oi) 











FIG. 1. Color online: Comparison of energy as obtained for different 
axial stripe phases at (a) /i/t = 1.00 and (h) h/t = 1.10. Ao is 
the strength of the modulation while the modulation vectors are q = 
(q, 0) with q = 2n'K jL. We have compared these also with diagonal 
stripes and 2D patterns, those data are not shown here. 


where Oi is the phase of and Oq is the phase at a fixed ref¬ 
erence site on the lattice, (c) the magnetization variable rrii = 
— n^), and (d) particle number Ui = {ui-^ + n^). These 
explicitly highlight the modulated nature in the FFLO state 
and spatial fluctuations with increasing temperature, (iv) The 
momentum occupation number ((nka)) that carries the sig¬ 
nature of population imbalance and finite Q pairing. Finally, 
(v) The fermionic density of states (DOSj^. 


III. RESULTS: THE GROUND STATE 
A. Variational phase diagram 

We computed the energy of the SAF Hamiltonian for three 
families of trial configurations {A^}. These correspond to 
axial stripes, checkerboard modulation, and diagonal stripes, 
with different wave vectors q, as detailed before. For a 24 x 24 



FIG. 2. Color online: Field dependence of magnetization (m) and 
axial stripe wave vector {q) at /x = — 0.2t for varying h as obtained 
from the variational calculation. On our system size we observe q ~ 
3m, while a more elaborate analysis suggests q = 7rm. 


FIG. 3. Color online: Variational ground state in terms of (a) varying 
— h and (b) varying n — m, showing the different LO phases. The 
results are on a 24 x 24 lattice. The allowed q values are discrete, 
titt/L, and the transitions between LO phases have small density 
discontinuities on this lattice size. The variation in q is expected to 
be continuous with m as L ^ oo, and the density discontinuities 
would vanish. The ‘operating point’ for which most of our thermal 
data is shown in this paper is marked by a circle in both panels. 


system we find that the axial stripes have the lowest energy. 
We have checked our results at a larger system size of 48 x 48 
and have found that the axial stripe phases still turn out to be 
energetically favorable compared to the other variational can¬ 
didates. A similar calculation by Chiesa et observed that 
Situ = At diagonal and axial stripes make up the ground state 
in the relevant n — m window. The possible reason for the 
differences between our results and those presented by Chiesa 
et al. are discussed at the end of the paper. 

Fig. I shows illustrative data on the energy of variational ax¬ 
ial stripes, A^ = Aocos{qxi), with respect to the amplitude 
Ao for a given q. The full variational comparison involves di¬ 
agonal stripes and checkerboard patterns as well, we have not 
shown them for clarity. Panel (a) shows that the q = 0 USF 
state is clearly unfavored and the global minimum is obtained 
for q = {tt/S, 0}. Panel (b), which is at a higher magnetic 
field, prefers a larger wave vector, q = {57r/12, 0}. 

The energy was calculated deliberately for finite size sys¬ 
tems, L = 24 in this case, since the Monte Carlo is done on 
that size. The absolute minimum, in the space of q and Aq, 
defines the mean field ground state for a given fi and h. The 
local minima are metastable within the mean field scheme. 

Given the modulated nature of the candidate states there is 
a size dependence to some of the features. By varying L from 
16 to 60 we found that (i) the minimum ^(/x, h) and Ao(/i, h) 
are reasonably size independent for L > 20, but (ii) the de¬ 
tailed ‘energy landscape’, as in Fig.I, does depend on system 
size, and for L > 48 the local minima vanishes, leaving us 
only with the global minimum. 

The overall reliability of our phase diagram, even on L = 
24, is borne out by our observation q ^ 7rm, Fig.2, which 
was obtained elsewhere on an infinite system using a more 
elaborate mean field decompositioiP^. We show m and q for 
varying h in Fig.2. The regime of validity of the relation 
q rsj Trm has been numerically establishecP^ and its basis has 
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FIG. 4. Color online: Spin resolved and total density of states for two different modulation wave vectors calculated through the BdG (red 
dotted line) and the Green’s function (black solid line) formalism. The dotted curves are shifted along the y-axis for the sake of clarity, (a) 
and (c) correspond to the spin up DOS for q = {tt/S, 0} and {vr/d, 0} respectively, while (b) and (d) correspond to the total DOS for the 
same modulation wave vectors, (e) and (f) shows the spin up spectral function A^^(k, uj) calculated for the {0,0} to {tt, tt} momentum scan 
through the Green’s function and BdG formalism, respectively, for q = {tt/S, 0}. 


been analytically discussed^ in the literature, in the context 
of half-filling repulsive Hubbard model. The repulsive Hub¬ 
bard model is related to the attractive Hubbard model via a 
particle-hole transformation. 

The q r-j Trm relation is based on two inputs: (i) the knowl¬ 
edge that the repulsive half-filling problem has (tt, tt) order, 
and (ii) the assumption that slightly away from half-filling the 
Fermi surfaces still retain the tilted square shape that they had 
at half filling. The second assumption is not obviously valid 
at large magnetization in the attractive problem, even if the 
mean density is n = 1, and makes the approximation unreli¬ 
able. Nevertheless it establishes a simple benchmark at small 
polarization. 

Fig.3 shows the ji — h and n — m phase diagrams inferred 
from the variational calculation. We have demarcated the re¬ 
gions pertaining to the different wave vectors. There is no 
noticeable window of phase separation in the n — m plot. The 
number density shows small discontinuity between the phases 
with different modulation wave vectors, as a consequence of 
the finite size of the lattice. The transitions are essentially con¬ 
tinuous and we expect the q variation to become continuous as 
L —cxD. 


B. Density of states 

Unlike the homogeneous superfiuid the modulated state 
gives rise to additional features in the energy spectrunJ^®!. 
Fig. 4 shows the spin up and total DOS as calculated through 
BdG diagonalization in comparison to those obtained through 
Green’s function formalism. These will serve as a reference 
when examining the thermal effects. The DOS deviates con¬ 
siderably from the homogeneous superfiuid case, and when 
finite size effects are el imina ted one expects a DOS with even 
more intricate featureJ^^® as shown through the Green’s 
function results. The Green’s function shows how the effect 
of scattering from to — k + and — k — in presence 

of the order parameter modulation. 

For any finite Q, the DOS deviates considerably from that 
of its BCS (Q = 0) counterpart. There is finite DOS at cc = 0 
and there are visible new van Hove singularities. These arise 
from the k regions where the dispersion Ea(k) satisfies the 
condition dEa{l<.)/dk = 0, where, a correspond to the mul¬ 
tiple branches that arises in the dispersion of the LO state. 

In Fig.4f we show the spin up spectral map for a (0, 0) ^ 
(tt, tt) momentum scan calculated through the BdG formal¬ 
ism, and compare the same with the result from the Green’s 













































































function approach, 4e. The dispersion shows three distinct 
branches separated by ‘gaps’ These regions of suppressed 
weight (not isotropic on the Brillouin zone) lead to the princi¬ 
pal depressions in the DOS shown in Fig.4a and 4c. 


IV. RESULTS: THERMAL BEHAVIOR 
A. Phase diagram 

Our first paper presented the overall h — T phase diagram at 
U = At. We observed that the unpolarised superfiuid ground 
state, with q = 0, undergoes a first order transition to a q / 0 
LO state at some h = hM The modulation wave vector 
of this state, and its magnetization, grows with field till, at 
h = hc 2 , the order is lost to a partially polarized Fermi liquid 
through a second order transitioJ^. For the U and fi that we 
have chosen, hd ^ 0.85f while hc 2 ^ 1.25t. hd and hc 2 
define the field boundary in the phase diagram in Fig.5. 

Fig.5(a)-(b) presents the — T phase diagram determined 
from our Monte Carlo calculation. The left panel shows the 
LO window estimated from the cluster based MC, with the 
superposed dotted line indicating the estimated from a 
small size exact diagonalization (ED) based MC. We have 
performed this check since the energy of small q LO states 
is sensitive to system size and cluster based approximations 
can give an erroneous estimate of Tc. The consistency with 
ED based MC suggests that our Tc estimates are reasonable. 

The typical Tc scales are ^ 0.0It. The wide ‘LO fluctu¬ 
ation’ window above Tc is characterized by the presence of 
finite q signature in the pairing structure factor, 5'A(q), and 
extends to T ^ O.lt near hd. We will show the 5'A(q) and 
r(q) further on and note that in the ‘balanced’ superfiuid such 
fiuctuation would be centered at q = 0. 

Fig.5.(b) compares the MC inferred Tc to the mean field es¬ 
timate. The MF Tc differs from the MC estimate by a factor 
of ^ 6 at hd, and a factor > 10 in the middle of the LO win¬ 
dow. The temperature axis in Fig.5.(b) is logarithmic. The 
lower two panels in Fig.5 show the MC based order parame¬ 
ter (on the left) and the MET based order parameter (on the 
right). They both involve strong first order transitions, but the 
Tc scale for MET is far higher as we have already noted. 


B. Structure factor 

In Fig.6 we show the thermal evolution of the structure fac¬ 
tor associated with pairing. We have explored both heating 
(from the variational state) and cooling (from an uncorrelated 
high T state). On heating from the variational ground state 
the axial stripes disorder (as would be more obvious from 
the spatial patterns later) and the structure factor peaks at 
(0, TQ) weaken and the system jumps discontinuously to the 
unordered but weakly correlated phase. The disordered phase 
has a fourfold symmetric feature in the structure factor for T 
just above Tc and this distorts into a more diagonal pattern by 
the time T ^ lOTc and is replaced by the usual broad feature 
around q = (0, 0) by T > O.lt. The highest T plot shows this 


feature. The fermionic correlation r(q) roughly follows the 
behavior of 5'A(q)- 

On cooling the pairing structure factor retraces the heat¬ 
ing features down to T ^ Tc but then transits to a some¬ 
what different, and poorly ordered, modulated state of the 
form Ai ^ e'^^^cos{Qy) at low T. This state is energeti¬ 
cally higher than the variational ground state, ^ cos{Qy). 
We discuss about this observation below. 

In the presence of metastable states the result of changing 
parameters like field or temperature can be path dependent. 
This is particularly true of Monte Carlo calculations where 
the system is supposed to explore the whole energy landscape 
and, for extended time, may remain trapped in local minima. 
As a consequence the ’heating’ and ’cooling’ cycles might 
lead to different low temperature states as shown in Fig.6. 
This is a computational observation, within the limits of sys¬ 
tem size, timescales and update method used by us. 

The major difference, with respect to this in the real system, 
would arise from two sources (i) the much larger size leads to 
a different pattern of metastability, and (ii) timescales: Monte 
Carlo accesses about 10^ sweeps per temperature, and the sys¬ 
tem can remain trapped on that timescale, while real life sys¬ 
tems involve timescales ^ 10“^s, i. e. ^ 10^ microscopic 
moves per second and may be able to escape metastable states 
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FIG. 5. Color online: h — T phase diagram and order parameter, 
(a). MC based h — T phase diagram at /i = —0.2 obtained on a 24 x 
24 lattice through heating cycle. It indicates the LO ordered and LO 
fluctuation regimes (see text). The dashed line corresponds to the Tc 
inferred from a 10 x 10 ED based MC calculations, (h). h — T phase 
diagram comparing the MC and MFT based Tc estimates. The T 
axis is logarithmic, (c) The temperature dependence of the ordering 
peak in Sa (q) at different h. (d) The temperature dependence of the 
MF order parameter. The temperature range in (c) and (d) differ by 
more than a factor of 10. 
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FIG. 6. Color online: Thermal evolution of the structure factors associated with pairing, 5'A(q) and r(q). The x and y axes in the panels 
represent Qx and Qy, both ranging from —it to tt. The upper set of panels is for heating from the variational LO state and the lower set for 
cooling from a random state. The structure factors display a q / (0,0) centered feature well above Tc, upto T ~ 0.05t. The cooled system, 
at the lowest T is unable to organize itself into a perfectly ordered axial phase. 


on laboratory timescales. 

The associated magnetic structure factor, not shown here, 
shows a prominent {0,0} peak, corresponding to bulk mag¬ 
netization, and a feature at finite q arising from ripples in the 
rrii, induced by the modulation in A^. This ‘antiferromag¬ 
netic’ peak is usually the only direct signature of a LO state. 
Neutron scattering experiments carried out on Pauli limiting 
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FIG. 7. Color online: Temperature dependence of (a) magnetization 
for varying h, and (b) the AF peak in the magnetic structure factor at 
h = 1.0 


superconductors do show such a featurJ^. Our finite size cal¬ 
culation makes it difficult to establish the lineshape of the AF 
peak, and anyway above Tc this peak is too faint to be visible 
compared to the bulk magnetization. 

Fig.7(a) shows the T dependence of the net magnetization, 
m = {ui^ — riii), at different magnetic fields. At lower h 
the system shows a first order transition at T = and for 
h ^ l.lt the m(T) is smooth - consistent with a very low 
ordering scale (our hc 2 estimate is ^ 1.2t). Fig.7(b) shows 
the temperature dependence of the antiferromagnetic peak in 
S'^(q). The discrete momenta on the 24 x 24 lattice makes it 
difficult to extract a meaningful lineshape for the AF peak. 


C. Spatial behavior 

Apart from the structure factor, the FFLO state can be char¬ 
acterized through real space signatures viz. pairing field am¬ 
plitude, phase correlation, magnetization and number density 
distribution. We present this in Fig.8. The figure clearly shows 
the axial stripe character of the paired state with the first two, 
low temperature, rows revealing that the modulation in | | 

and rrii have the same period while the phase variable has dou¬ 
ble the wavelength. The ‘nodes’ in |A^| correspond roughly 
to the peaks in m^. The local density rii remains within 1% of 
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FIG. 8. Color online: Spatial maps characterizing the thermal evolution of the LO state at /i = l.Ot, through (a) pairing field magnitude 
(| Ai |), (b) phase correlation (cos(^o — is a reference site), (c) magnetization (rrii = — n^) and (d) number density (rii = + n^). 
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FIG. 9. Color online: Thermal evolution of the momentum occupation number na(k.) at h = l.Ot. 


0.94. The third row is for T > Tc and there are no obvious 
signatures of even short range correlation in the amplitude and 
phase variables, while rui is essentially homogeneous. Nev¬ 
ertheless, as S/\{q) reveals, weak finite q correlations survive 
to a fairly high temperature. 


D. Momentum distribution 

We compute the spin resolved fermion momentum distribu¬ 
tion and show the result for the heating cycle dXh = tm Fig. 9. 
The symmetry is clearly two fold due to the axial modulation. 


In the FFLO state the pairing is between k, t and —k + Q, 
where the k’s are supposed to be on the up Fermi surface (FS) 
and the — k + Q on the down FS. The Q inferred from the 
structure factor and the spatial maps is ±^7r/3 for the axial 
LO state. 

To check how well the paired momenta lie on the respec¬ 
tive FS we varied k on the up spin FS and located the corre¬ 
sponding —k ± ^tt/S on the down spin FS. For k in the upper 
branches of the up spin FS the connection generates 

points that are reasonably close to the smaller down spin FS. 
The correspondence is best for points along the diagonal and 
worsens as we move k along the FS to the zone boundary. The 
































































FIG. 10. Color online: Density of states and the distribution function for | A|. (a). Temperature dependence of the up spin DOS T) at 

h — t, (b) Temperature dependence of the total DOS N(uj^ T) at /i = t, and (c) P(| A|, T) at /i = t. The thermal evolution wipes off the 
delicate features of the low temperature DOS, associated with Andreev reflection. 


complementary behavior holds for the lower branches of the 
FS when we consider the — ^tt/S momentum transfer. 

For k neighbourhoods where k,t fails to find a partner 
as — k + Q,| the gap is smaller, and a seemingly nodal 
gap structure, arising due to the LO modulation, can be 
seen via the spectral weight distribution. We have discussed 
about these spectral features of the modulated LO phase, 
elsewherd^. 

Now the thermal evolution. For T < Tc ^ 0.005t the 
distributions show only a two fold symmetry and of course 
differing sizes due to the magnetization in the LO state. For 
T > Tc and up to ^ 2Tc the distributions are almost four¬ 
fold symmetric but have some weak modulation. By the time 
T ^ 0.05t ^ lOTc the distributions only refiect the popula¬ 
tion imbalance of the PPFL, with no trace of modulation. 


E. Density of states 

Although most of the ‘sightings’ of a FFLO state are in the 
solid state, there seem to be no detailed data confirming its 
unusual spectral features. In continuum atomic gases promi- 
nen t pseu dogap features have been observe d in the balanced 
cas jilEl] and analyzed in detail theoreticall} )^^^! For the im¬ 
balanced gases also both experimental and theoretical studies 
suggest pseudogap features at large imbalanceJ^33®, even 
though the detailed angle resolved measurements have not 
been made. Finite T results on the FFLO state on a lattice 
seem to be rare. We discuss the ‘pseudogap’ features in this 
regime based on our density of states observations, below. 

We focus on the following: (i) The low temperature spec¬ 
trum that is characteristic of the modulated pair state, with its 
multiple peaks and dips, (ii) The thermal evolution of the ideal 
T = 0 spectrum in the T < window, and (iii) The signa¬ 
tures of large U, and possibly short range pairing correlations, 
in the T > window. Note that in the Hubbard model, i.e, 
with contact interactions, the effectiveness of Fermi-Fermi in¬ 
teractions weaken with increasing population imbalance - so 
the survival of a high T pseudogap at large population imbal¬ 


ance is not obviouJ^. 

Fig.I0(a)-(b) shows the spin up and spin summed fermion 
density of states forh = t corresponding to Q = |Q, 7r/3|. 
The up spin spectrum in the ground state has a depressioiF^ at 
uj ^ —h, i.e, the LO ground state is ‘pseudogapped’ (PG), 
due to its peculiar band structure, unlike the gapped ‘BCS 
state’E^m Increasing T quickly weakens these features, and 
the strong first order transition at Tc leads essentially to a gap¬ 
less state. Tc therefore is also the PG to ‘ungapped’ transition 
point for the spin resolved DOS. The gapless phase survives to 
a scale beyond which the system again shows a weak 

pseudogap. This ‘re-entrant’ feature is due to the thermally 
induced large amplit ude fl uctuations in this strongly interact¬ 
ing, U = At, problenPSSI. 

The observed spectral features can be analyzed based on 
the distribution P(|A|). In Fig. 10(c) we have shown the dis¬ 
tribution of the pairing field at temperatures corresponding to 
those of the DOS. For T < Tc P(|A|) is multi-peaked due to 
the amplitude modulated LO state. Increasing T beyond Tc 
makes the distribution broad, single peaked, and reduces the 
mean magnitude of A^. The weak pairing field leaves the DOS 
essentially featureless. For T > 0.05f the behavior of the dis¬ 
tribution changes significantly. The mean of P(|A|) shifts to 
significantly high values with increasing T, although the dis¬ 
tribution remains very broad. This, we believe, is similar to 
the PG effect seen at large imbalance in the unitary Fermi gas. 


V. DISCUSSION 

We have discussed the strength and limitations of the static 
auxiliary field scheme in Paper There are size limitations 
that are specific to the FFLO phase which we touch upon here. 
We also discuss experimental situations where our results can 
shed some light. 














A. Methodological issues 

1. Single channel decomposition 

We have decomposed the Hubbard interaction only in the 
pairing channel. At T = 0 our auxiliary field approach, using 
the field A^, reduces to a ’pairing channel’ mean field theory, 
the standard BdG scheme. Mean field theory can be imple¬ 
mented in more complex formats, by decomposition in the 
pairing, density, and spin channels simultaneously. This is 
the Hartree-Fock BdG (HFBdG) approach and is obviously 
more general than simply BdG. Such calculations yield diag¬ 
onal stripes over a part of the phase diagranJ^I while we obtain 
only axial stripes. 

Unfortunately, the multichannel mean field theories are not 
good templates for including fiuctuations beyond the Gaussian 
level. Since the strength of our method is really in calculat¬ 
ing the thermal properties we opted to use the single chan¬ 
nel Hubbard-Stratonovich so that large amplitude fiuctuations 
can be systematically included, and the normal state well cap¬ 
tured. 


2. Quantum fluctuations in the FFLO state 

We use a static auxiliary field technique where the temporal 
(quantum) fiuctuations of the pairing field are neglected, while 
spatial fiuctuations are completely retained. While this could 
be a poor approximation in the continuum, on a lattice it is 
reasonable as we argue below. 

The low energy fiuctuations in the continuum FFLO state, 
i.e, in free space, arise from three sources (i) the ‘phase sym¬ 
metry’ of the f/(l) order parameter, (ii) the translational sym¬ 
metry breaking, and (iii) the rotational symmetry breakin^^. 
As a result, in two dimensions, long range order cannot be sus¬ 
tained even at T = 0 and mean field theory (which predicts 
such order) is invalid. 

On a lattice the phase field still has low energy excitations 
of the ‘XY’ type, but the translational and rotational modes 
would be gapped out due to the spatial symmetry already bro¬ 
ken by the underlying latticJ^. Models with XY symmetry 
have long range order in the 2D ground state and a KT transi¬ 
tion at finite T. Therefore, the issue of fiuctuations reduces to 
checking how well the U{1) superfiuid Tc is captured by our 
model, in 2D, vis-a-vis full QMC. A comparison between the 
QMC results and that obtained by the present numerical tech¬ 
nique has been made in ref.® However, for a spin imbalanced 
system such QMC results do not seem to exist. Nevertheless, 
the argument above and the Tc comparison for the balanced 
system shows that, for the fiuctuations that are relevant, our 
method does quite well. 


3. Finite size effects in the FFLO phase:- 

The variational ground state calculation for the FFLO state 
is affected by the finite size of the lattice. The number of q 


values (modulation wave vector of the LO state) is determined 
by the size of the lattice. 

A more serious size limitation arises for the finite tempera¬ 
ture Monte Carlo simulation owing to the cluster update tech¬ 
nique that we us^. Even though the cluster update technique 
allows us to access large system sizes (^ 30 x 30) within rea¬ 
sonable computation time, it introduces another length scale 
in the problem. Consequently, within our approach there are 
two sources of error, (i) The finite wavelength. A, of the FFLO 
modulations require the linear dimension of the system, L to 
be much greater than A. (ii) If the energy cost of MC update is 
not computed via exact diagonalization of the L x L system, 
but on a Lc X Lc cluster, one needs Lc ^ X. These are dif¬ 
ficult constraints to satisfy when trying to access large A, i.e, 
small Q states. As a result most of our detailed thermal results 
are on relatively large Q states. 


B. Connection to experiments 

As discussed in the introduction the FFLO state has so 
far been realized experimentally only in the solid state, with 
some evidence in heavy fermionJ^®® organic^®® and iron 
pnictide^I^®^. However, most of the superconductors are non 
5-wave, and are at considerably weaker coupling than we 
have considered. None of them are close to their BCS-BEC 
crossover coupling, so a detailed comparison of our thermal 
results with them is inappropriate. 

Turning to cold atoms, the key signature of a LO state is 
the spatial modulation of the pairing field and polarization. 
Moreover, the LO state has finite polarization down to zero 
temperature, unlike the breached pair state where the polar¬ 
ization vanishes as T 0. While there are so far no ex¬ 
perimental signatures of real space polarization modulation, 
the spin resolved density profiles in quasi one dimensional 
traps indicate finite polarization down to the lowest accessi¬ 
ble temperaturJ^. The spectral features of this state has not 
been measured. 


VI. CONCLUSION 

We have studied the attractive Hubbard model in the cou¬ 
pling regime of BCS-BEC crossover and large population im¬ 
balance. We established the thermal properties of the EELO 
phase that occurs in this regime. The thermal transition is 
strongly first order with a Tc much below the mean field 
prediction. Einite momentum pairing correlations neverthe¬ 
less persist to a scale that is ^ lOTc. We have studied the 
fermionic density of states and discover thermally induced 
closure of the characteristic EELO dip at T ^ Tc but the reap¬ 
pearance of a correlation induced pseudogap at higher tem¬ 
perature. The Tc estimate and the detailed physical indica¬ 
tors should help in pinpointing EELO states in strong coupling 
cold atomic optical lattices. 
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